Mon. Not. R. Astron. Soc. 000, ITlfl7l f2222) Printed 13 January 2010 



(MN MfeK style file v2.2) 



On the use of structure functions to study blazar 
variability: caveats and problems 



D. Emmanoulopoulos 1 *, I. M. M c Hardy 1 and P. Uttley 1 

1 School of Physics and Astronomy, University of Southampton, Southampton SOU 1BJ 



Accepted 2009 Month #. Received 2009 Month #; in original form 2009 Month # 



ABSTRACT 

The extensive use of the structure function (SF) in the field of blazar variability across 
the electromagnetic spectrum suggests that characteristics time-scales are embedded 
in the light curves of these objects. We argue that for blazar variability studies, the 
SF results are sometimes erroneously interpreted leading to misconceptions about the 
actual source properties. Based on extensive simulations we caution that spurious 
breaks will appear in the SFs of almost all light curves, even though these light curves 
may contain no intrinsic characteristic time-scales, i.e. having a featureless underlying 
power-spectral-density (PSD). We show that the time-scales of the spurious SF-breaks 
depend mainly on the length of the artificial data set and also on the character of the 
variability i.e. the shape of the PSD. 

The SF is often invoked in the framework of shot-noise models to determine the 
temporal properties of individual shots. We caution that although the SF may be fitted 
to infer the shot parameters, the resultant shot-noise model is usually inconsistent with 
the observed PSD. As any model should fit the data in both the time and the frequency 
domain the shot-noise model, in these particular cases, can not be valid. 

Moreover, we show that the lack of statistical independence between adjacent SF 
points, in the standard SF formulation, means that it is not possible to perform robust 
statistical model fitting following the commonly used least-squares fitting methodol- 
ogy. The latter yields uncertainties in the fitting parameters (i.e. slopes, breaks) that 
are far too small with respect to their true statistical scatter. Finally, it is also com- 
monly thought that SFs are immune to the sampling problems, such as data gaps, 
which affects the estimators of the PSDs. However we show that SFs are also troubled 
by gaps which can induce artefacts. 

Key words: methods: statistical - methods: numerical - methods: data analysis - 
galaxies: individual: Mrk501 - galaxies: active - BL Lacertae objects: general 



1 INTRODUCTION 

The flux variation of Active Galactic Nuclei (AGN) is a 
phenomenon that can provide us with significant informa- 
tion about the physical properties of these sources. The 
most variable AGN are the blazars which exhibit dramatic 
flux variations across the electromagnetic spectrum. The 
fastest variations are observed in the X-ray and 7~ray bands 
on time-scales of hours or even minutes dCatanese et al.l 
ll997l : lMaraschi et al.l l l999l: lAharonian et al.ll2003l, 12005a 



whereas in the optical (iTosti et aLlll998l: [Stalin et al.ll200 
and radio wavebands llRomero et al.lll997l ; lAller et al.ll 19991 : 
iTerasranta et al1l2005l ') the temporal variations are seen on 
time-scales of days to weeks, up to years. 

In the last fifty years several time-series analysis meth- 
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ods have been developed to study the var iability properties 
of as tronomical sources (see for a review ISubba Rao et all 
Il997h . Linear and nonlinear analysis methods in the time 
and frequency domains, adjusted to the needs of astronom- 
ical data sets i.e. taking into account measurement errors 
and data-gaps, can describe adequately the time behaviour 
of AGN. Nevertheless, several authors, as we are going 
to discuss further, are convinced that more conventional, 
intuitive-tools, based on "running variance" computations 
are able to give robust results with respect to the underly- 
ing variability properties of the observed source. Although 
these methods have value in some cases, they can often give 
misleading results as we shall show in this paper. 

One of the most extensively used tools in the field of 
blazar variability i s the structure function (SF) (see e.g. 
iHughes et ailll992f) which measures the mean value of the 
flux- variance for measurements, x(t), that are separated by 



2 Emmanoulopoulos et al. 



a given time interval, r, where SF(t) = ([x(i) — x(t + t)] 2 ). 
The SF is commonly characterized in terms of its slope /3, 
where SF(t) oc t 13 . 

Historical ly the SF has been u se d in the study o f turbu- 
lent plasmas i|Kolmogorovlll941a! lbT; lYu et all 120031 ). It was 
introduced to astrophysics by radio astrono mers studying 
slow scintillations in the interstellar medium l|Rickett et all 
1 19841 ) following the methodo logies o f lProkhorov et all (|l975l ) 
and IColes fc Frehlich] (|l982l ) on the subjects of laser prop- 
agation in turbulent media and atmospheric scintillation 
respectively. The first systematic description of the SF 
methodology, adju sted to the needs of a stronomical data 
sets, was made by Simonetti et all (|l985l ) using as a refer- 
ence the work of iRutmanl ( 19781) from the fi eld of electrical 
and electronic engineering. Simonetti et al.l used the SF to 
demonstrate that the time-series of flat- and steep-spectrum 
radio sources differ qua litatively. During the same period, 
Cordes fc Downs! (Il985l) esti mate the SFs of 21 pulsars and 



Hiellming fc Naravanl ( | 19861 ) derived the first quantitive SF 



results for the compa c t galactic radio s o urce 1741-038. Then, 
iFiedler et al. | (| 19871 ); iHeeschen et"aH (|l987l ) also used the 
SF to study the variability properties of large extragalac- 
tic radio samples. Subsequently, the SF has been employed 
for the study o f the timing pro p erties o f objects in higher 
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energy bands. Bregman et al.l (fl988|); Neugebauer et al 



Quirrenbach et al. ( 1992) ; iHufnagel fc Bregman 
Smith et alj i 19931 ) ; iLainela fc Valtaoial (l993l) 



Heidt fc Wagnerl (Il996[) w hilst iLewis fc Irwinl (|l996l ) and 



Bland ibrd fc Kundid ( 19971 ) applied the SF-analysis method 
to derive masse s for microlenses. In the X-ray band, 
iBrinkmann et al. I (|2000l ) applied the SF methodology to 
the ROSAT observations of PKS 2155-304 and they at- 
tributed the derived shot-noise variability to either the ac- 
celerating inner jet model or to accelerating particles at 
shocks travelling down the relativistic jet. Thereafter many 
X-ray variabilit y many X-ray observations have employed 
SF analysis (e.g. Gliozzi et al.l 2001 Brinkmann et alj|200ll ; 



llvomoto fc Makishimall200ll ; IZhang et al.l l2002) in order to 
derive some short of characteristic time-scale. 

SFs have been particularly commonly applied to obser- 
vations of blazars. The ASCA "long -look" observations of 
Mrk501, Mrk421, and PKS 2155-304 l|Takahashi et al.ll2000l: 
iTanihata et al.ll200ll . l2003l ) are, even now, some of the most 
uniformly sampled light curves in the X-ray regime. The 
main outcome of the above SF analysis is an apparent char- 
acteristic time-scale of ~1 day which appears to be in accor- 
dance with the duration of the shots ex isting in the inter- 
nal shock particle acceleration scenario (jSikora et al.ll200ll ; 
ISpada et al.ll200ll ). 

In this work we show through extensive simulations that 
much of the existing literature on blazar SFs tends to misin- 
terpret observed SF characteristics such as breaks and slopes 
as being real or physically meaningful, when often these are 
either artefacts intrinsic to SFs, or subject to much greater 
statistical variation than inferred from the commonly-used 
fitting procedures. In Section[2]we present the method that 
we are going to use in order to create artificial light curves 
lacking any sort of characteristic time-scales and in Section 
[3] we study the behaviour of the SF derived from the thor- 
ough ly studied ASCA data set of Mrk501 (|Tanihata et all 
l200ll . hereafter TAN01). For the same source we employ 
the long- look light curve of All-Sky Monitor (ASM), onboard 



Rossi X-ray Timing Explorer (RXTE), in order to study the 
effects of the data length on the position of the SF-break. 
Then, in Section U we study the timing properties in the 
Fourier domain of the shot-noise model that TAN01 have 
used in order to study the SF properties. In Section [5] we 
test the statistical robustness of the most commonly used 
fitting procedures which are employed in order to derive as- 
trophysically interesting quantities from the SF. Finally, in 
Section [S] we study the sensitivity of SF to the presence of 
data-gaps. 



2 LIGHT CURVE SIMULATIONS AND 
POWER SPECTRAL DENSITY 

We simulate stationary light cu r ves b ased on the proce- 
dure used by iTimmer fc Koenig] (|l995l ) which uses as an 
input the Power Spectral Density (PSD) function of the 
observed light curve and gives an ensemble of stochastic 
time-series produced from the same underlying PSD. This 
is the most robust method of mimicking the properties of 
a light curve since it takes into account correctly the in- 
trinsic scatter in the power at a given frequency, described 
by a xl, chi-squared dis tribution with two degrees of free- 
dom (e.g. |Priestlevl ll981). by randomizing both phases and 
amplitudes. 

Throughout this paper we use the standard nomen- 
clature, i.e. we refer to the periodogram as the modulus 
squared of the discrete Fourier transform of the partic- 
ular data set under consideration whilst the PSD repre- 
sents the mean of the true underlying dist ribution of vari- 
ability power as a function of frequency (|Priestlevl Il98ll ; 
IVaughan et~ai1l2003l ). Since the periodogram is highly scat- 
tered aro und the PSD hav ing a standard deviation of 100 
per cent (jPress et al.lll992D. we use the binned logarithmi c 
periodogram (|Papadakis fc Lawrencell 19931 ; IVaughan|[2005l ) . 
If there are more than 20 samples in each geometric-mean 
frequency bin then their distribution will be Gaussian, thus 
producing a statistically sensible standard deviation. 

Concerning the length of the simulated light curves 
some special caution should be taken with respect to the 
effect of red-noise leak i.e. the transfer of power from low to 
high frequencies by the lobes of the window function which 
can produce features such as slowly incr easing or failin g 
trends across the generated light curve (e.g. |Priestievlll98ll ). 
To take this effect into account, we extend the PSD to very 
low frequencies (i.e. long time-scales), in order to generate 
artificial light curves 50 times longer than the observed one. 
We then truncate the simulated data trains to the desired 
length by selecting a random segment having equal length 
to the real light curve under study. 

In the case of blazar variability the underlying PSD can 
be very-well described by simple or broken-power-laws with 
indices 1 < a < 2.5 flattening at long time-scales depict- 
ing the physical fact that the variability amplitude can not 
increase forever (stationarity) . For these kinds of physically 
stationary processes the periodogram can correctly describe 
the underlying PSD even for steeper slopes of -4 (or even 
-5). 

Apart from red-noise leak, spectral representations are 
also affected by aliasing effects. Frequency components out- 
side the frequency range covered by the data set, are 
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aliased into that range by the very act of discrete sam- 
pling (|Press et al.l |l992). Finally, irregular sampling pro- 
duces a window fun ction which smears the spectral esti- 
mates (|ScargleHl989Tl . causing deviations from the true un- 
derlying PSD. Over the years a number of methods have 
been derived to overcome these problems. In particular, 
we can estimate t he underlying PSD by model fitting (e.g. 
iDone et all Il992l ; lUttlev et all l2002h . This process is now 
commonly appli ed to observations of Seyfert galaxy X-ray 
variability (e.g. iMarkowitz et al.l|2003l ; lUttlev fc McHardvl 
l2005l ; lMcHardv et al.ll2004l2005l , [20061 , |2007f l and produces 
reliable results. Similar simulation-based treatments of the 
SF results are missing from the blazar literature. 



3 CAVEATS REGARDING THE STRUCTURE 
FUNCTION BREAK-TIME-SCALES 

3.1 The ASCA data set of Mrk501 

To study the behaviour of the SF, we examine the ASCA 
data set of Mrk501 (TAN01) in 2-10 keV, with a sample 
interval of 5678.3 sec. This is one of the three data sets (the 
other two are Mrk421 and PKS 2155-304) in which TAN01, 
based on the SF analysis method (Fig. [1] left panel), find 
characteristic time-scales and suggest that these are a sig- 
nature of the minimum time-scales of individual shots. This 
is one of the longest and most continuous data sets ever ob- 
tained for a blazar in the X-rays, covering a time span of 10 
days. 

Initially, we estimate the binned logarithmic peri- 
odogram of the data set which can be very well fitted by 
a simple featureless power-law with index a = —1.80 ± 0.09 



of the break. Even if a clear plateau is not formed after 
t > Tbr, several authors consider the first "hump" as a sig- 
nature of a genuine t e mporal property of t he source (e.g. 
Takahashi et ail l20od iGliozzi et al.l l200ll; iKataoka et all 



for 6 degrees of freedom with a null hypothesis 



probability of 0.70). There is no evidence for a break in the 
PSD. 

We next produce 2000 artificial light curves having the 
same power-law PSD of index a — —1.80 and the same 
length as the studied ASCA data set. For each artificial light 
curve we estimate the SF by employing exactly the same 
functional form of the SF as TAN01 

SF{r) = E W (*M* + + T ) ~ /»] 2 (!) 

where N(t) — ~}2 w(i)w(i+r), r is the separation time of the 
observations, and w(i) is the weighting factor which is pro- 
portional to the data point f(i) and inversely proportional 
to its error 07(1). By an eye inspection, we note apparent 
breaks in almost all of the simulated SFs. 

In our simulations we take into account the effect of 
the measurement errors in the SF-estimates by represent- 
ing each measurement as a deviate drawn from a Gaussian 
distribution with mean and standard deviation equal to the 
measurement's value and error respectively (200 times for 
each artificial light curve). We have to note that for this 
particular ASCA data set the measurement errors do not 
play a major role in the SF-estimates. If we ignore these er- 
rors the values of the SF-breaks and SF-slopes change only 
by ~ 1 per cent and ~ 0.5 per cent respectively. 

Next, we determine the position of the SF-break (Fig. 
[T] right panel). To localize the break we produce an in- 
terpolated version of the SF and we find the abscissa Tt r 
of the first local maximum, which indicates the position 



200ll ; lAgudo et alj|2006t iFuhrmann et al.ll2008l ) 

Finally we note that after the first break point, the SF 
usually exhibits "wiggly-patterns" forming a plateau. These 
fluctuations indicate that the pairs {f(i + t), f(i)}r>r bT 
which are averaged within each time bin r > Tb r , are lin- 
early independent having zero linear correlation (Appendix 
[A")) . In order to study the frequency of occurrence of these 
SF structures, we register the abscissas of all the local ex- 
trema occurring for r > Tb r . 

In the left panel of Fig. [2] we present with the filled 
grey area the distribution of the SF-breaks coming from our 
simulated light curves that have the same length and the 
same featureless PSD as the original ASCA light curve of 
Mrk501. From an ensemble of 2000 simulated data sets, we 
get breaks from 1903 of them. To specify the position of 
the maximum in our histogram we fit the histogram entries 
with a Gaussian distribution. The results of the fit yield 
an amplitude of 279.55 samples, a mean value of 0.91 days 
and a standard deviation of 0.09 dayfl With the solid and 
dotted lines we present the distribution of the local extrema 
(maxima with solid line and minima with the dotted line) 
occurring after the first break. The occurrence times of the 
"wiggling" patterns are purely randomly distributed for r > 
Tbr following a uniform distribution. 

These simulations show us clearly that stochastic data 
sets having the same length as the ASCA data set of Mrk 501 
and the same featureless PSD, exhibit breaks in the SF 
around a day. Thus, the apparent break seen in the SF for 
Mrk 501 (Fig. [T] left panel) should not be associated with 
any sort of physically meaningful time-scale. 



3.2 Longer time-series: The ASM data set of 
Mrk 501 

Following Section 13.11 we extend our simulations to even 
longer time-scales to check how SF-break time-scales might 
be related to the length of a data set. Assuming that for 
long time-scales the variability properties of MRK 501 are 
well-represented by the same PSD describing the short-term 
variability (i.e. same slope and normalisation), we produce 
2000 artificial light curves, each being 2000 days long. For 
every light curve we calculate the SF and we specify the 
position of the SF-breaks. In total 1893 light curves exhibit 
SF-breaks and the distribution of their position can be well 
represented by a Gaussian distribution having a mean 399.35 
days and a standard deviation of 74.73 days (Fig. [2] right 
panel) . 

In order to compare our predicted break with the data 
we use the long-term RXTE-ASM light curve of Mrk 5010 
spanning from MJD 50100 to MJD 52100 and we estimate 



1 By fitting a log-normal distribution the results remain practi- 
cally the same giving an amplitude of 254.33 samples, fi = —0.09, 
and a = 0.10, yielding a mean value of 0.92 days and a standard 
deviation of 0.09 days. 

2 We have obtained the light curve from the ASM light curve site 
of MIT: |http : //xte .mit . edu/ASM_lc .html| 
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Figure 1. [Left panel] The SF (in logarithmic scale) for the ASCA data set of Mrk501 as shown in panel (a) of figure 4 in TAN01. The 
arrow indicates the break time-scale which is around a day. 

[Right panel] A randomly chosen SF (in logarithmic scale) from the ensemble of the 2000 artificial light curves which are produced from 
a featureless power-law PSD of index -1.80 with no characteristic time-scale, having the same length as the ASCA data set. There is a 
clear break around one day, similar to the one derived from the ASCA SF. 
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Figure 2. The distribution of the SF's breaks and wiggling features. 

[Left panel] The grey area depicts the distribution of the SF-breaks coming from an ensemble of 2000 artificial light curves which are 
produced from a featureless power-law PSD of index -1.80 with no characteristic time-scale, for the case of the ASCA data of Mrk501 
(the histogram bins have a length of 0.04 days). Based on a Gaussian fit the mean value and the standard deviation of the distribution 
are 0.92 and 0.09 days respectively. The solid and dashed lines represent the distribution of the local maxima and local minima for 
t > Tb r mapping the positions of the wiggling features. 

[Right panel] The distribution of the SF-brcaks for the case in which the simulated light curves, produced by the same PSD as above, 
extend to 2000 days (the histogram bins have a length of 10 days). It has a mean value of 399.35 days and a standard deviation of 74.73 
days. 
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its SF. The SF-break occurs around 402 days, something 
which is absolutely in accordance with the results from our 
simulated light curves of the same length which did not in- 
clude of any sort of characteristic time-scale. 

The aforementioned example reflects in a clear way that 
the SF deals only with the properties of the observed light 
curve ignoring the properties of the true underlying variabil- 
ity process. During astronomical observations we observe a 
source for a given time and from the observed data series 
we have to extract in a statistically robust way the proper- 
ties of the underlying variability process. If the result that 
we obtain is a strong function of the length of a given data 
set then this can lead to a serious misunderstanding of the 
variability properties of the source. 

Similar SF effects to those derived for the ASCA and 
ASM light curves of Mrk501 should be expected for all sim- 
ilar light curves. By producing an artificial light curve 6000 
time units (t.u.) long, from a featureless power-law PSD of 
slope -2, we estimate the normalized SF (NSF) (equation 
()A10f> ^ (i.e. the SF divided by two times the value of the 
variance of the data set, see Appendix [X| for the overall 
data train and for a 500 t.u. subset of it (Fig. [3]). In this 
case we use the NSF, instead of the classical SF, for an easy 
and direct comparison of the position of the break, since the 
formed plateaus in both SFs are distributed around 2. The 
top right panel of Fig. [4] shows the behaviour of the NSF for 
these two data sets, coming from the same underlying vari- 
ability process. The break in the NSF for the short data set 
occurs at around 40 days and for the overall light curve at 
1600 days respectively. In this case, the fact that the break 
occurs at two different time points for two different dura- 
tions of the data set which is fully described from the same 
and completely featureless PSD, indicates that the break 
does not reflect any physically interesting time property of 
the process itself. 

For comparison purposes we also estimate the binned 
logarithmic periodograms for the two light curves discussed 
above. As we can see from the bottom left panel of (Fig. [4]) 
both can very well reproduce the input PSD slope of -2. 

3.3 SF-breaks and data set length 

To quantify the relationship between the spurious SF- 
breaks, Tbr, and both the length and the PSD-slope of the 
data set, we perform a set of simulations. For each feature- 
less power-law PSD (with index a = — 1, a — —1.5, a = —2, 
a — —2.5 and a = —3) and for various data-set lengths (100, 
400, 700, 1000, 1300, 1600, 1900 and 2200 t.u.) we produce 
2000 artificial light curves. In total we perform 40 simu- 
lations, each one consisting of 2000 artificial light curves. 
Then, for each PSD-index and light curve length we localise 
the position of the Tbr, as described in Section [371] and after 
forming its distribution we estimate its mean and its stan- 
dard deviatiorjf] (Fig. [5j . 

For all the PSD-slopes we can readily see that as the 
length of the data set increases from t t.u. to t' t.u. (t' > t), 
the SF-break gradually shifts to larger values. This can be 
understood for the case of a random walk process (a = —2) 
where the mean distance after T steps is vT. The mean 

3 The distributions are very close to Gaussian (see e.g. Fig.[2]l. 



distance in our case is the mean variance which is measured 
from the SF, hence for t' > t the maximum variance and 
hence the position of the break is expected to shift towards 
larger time-scales. 

The important point here is that the measured variance 
of a data set (and hence the maximum variance depicted by 
the break) does not reflect the true underlying variability 
properties of the source. As the data set increases, larger 
variations produced from the same variability process having 
the same PSD enter the data sets and this is mapped in the 
SF as a break even in the case of featureless PSD which are 
lacking characteristic time-scales. 

In Fig. [S]we show the position of the SF-breaks, Tbr, 
induced by the different length of the data sets for various 
PSD-slopes. This diagram can be used only for continuously 
sampled light curves for which the underlying PSD is a pure 
power-law with the given indices. The existence of data- 
gaps in the data set can alter significantly this picture since 
sampling patters play a major role in the form of the SF 
and hence the position of the break (Section [6]). In the case 
that our data set has a PSD of a broken-power-law form, 
similar simulations should be performed in order to specify 
the position of the fake SF-break. Note that the physically 
interesting PSD-breaks will be mapped in the SF but in 
statistically unsound way with respect to their position and 
their uncertainties (Section 15. 3[) . 

Cons ider the case of the blazar PKS 2155-304 as pre- 
sented bv lZhang etafl (|2002r i. The 1997 (~100 points) and 
1999 (~200 points) data sets are the most continuous ones 
(figure 2 and figure 3 in Z hang et al.l 120021 ') and they have 
featureless PSD with slopes a = — 2 and a = — 3 respectively 
(Table 3 in lZhang et"alll2002l ). The SF of the corresponding 
light curves (figure 6, middle right pane l and bottom right 
panel respectively, in IZhang et al.l 120021 ) exhibit breaks of 
the order of Tb r ~ 12.8 ksec and Tbr ~ 61.7 ksec respectively 
(or 13 t.u. and 61.7 t.u respectively for SF-bins of 1 ksec). 
A SF-break for the 1997 data set is predicted (from Fig. [5} 
to lie at around 14.8 ± 1.4 t.u. and for the 1998 data set, at 
around 60 t.u. with a simple linear interpolated estimation. 
This outcome readily tells us that these breaks can originate 
from a variability process with no characteristic time-scale 
simply described by a featureless PSD of a power-law form. 



4 THE STRUCTURE FUNCTION AND THE 
BLAZAR SHOT-NOISE MODEL 

Particle acceleration in shocks (|Kirk et al.l Il998h is com- 
monly invoked to explain the observed temporal and spec- 
tral variability of blazars. There may be a variety of intrinsic 
timescales, i.e. acceleration, escape and light crossing time- 
scales, associated with a single emitting region. TAN01 sim- 
ulate light curves on the assumption of certain shot parame- 
ters and, from the similarity between the simulated and ob- 
served SFs, estimate some temporal shot properties. Here we 
caution against this approach as, in general, the light curves 
derived from such shots do not exhibit the same PSD as 
those observed from the blazars. A physically correct model 
should reproduce both the SF and the PSD. 

In TAN01 a number of simulations were carried out in 
order to associate the SF-breaks with the properties of the 
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Figure 3. [Left panel] The long-term 2-10 keV TiSM light curve of Mrk501, between 50100 MJD-52100 MJD, in bins of 15-days. The 
dotted line is linear interpolation intended to guide the eye. 

[Right panel] The SF of the ASM light curve of Mrk501, in bins of 15-days (in logarithmic scale). The SF-break occurs around 402 days. 
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Figure 4. [Left panel] The overall simulated light curve 6000 t.u. long produced from an underlying PSD having a power-law form with 
index -2. The inlay shows a short segment of the light curve 500 t.u. long. 

[Top right panel] The NSF of the data in the inlay (empty diamonds) and the overall data set (black points). The break of the SF occurs 
at 40 and 1600 t.u. respectively. 

[Bottom right panel] The binned logarithmic periodograms of the short (filled diamonds) and the overall data sets (grey points). Both 
of them match the underlying PSD, being a featureless power-law of index -2. 



observed flares. The light curve is regarded as a superposi- 
tion of triangular shot^\ with rise and decay time-scale t T 



and id respectively, occurring randomly at t p following a 



4 The last 40 years the commonly used shot-noise models consist file or only an expo nential decay profile (e.g. lLochner et al.l[l99ll : 

of pulses with either a symmetric exponential rise and decay pro- iBurderi et al.|[l997l) . 
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Figure 5. The position of the SF-break, Tb r (in logarithmic 
scale), as a function of the length of the data set assuming an 
underlying PSD of a power-law form of index a. Each data point 
represents the mean value and the standard deviation of the dis- 
tributions of the SF-breaks yielding from an ensemble of 2000 
simulations having a power-law PSD of index -1 (filled circles), 
-1.5 (filled squares), -2 (filled diamonds), -2.5 (filled up-pointing 
triangles) and -3 (filled down-pointing triangles). The dashed lines 
are linear interpolations intended to guide the eye. 



Poisson distribution, with an intensity Io at t p 

( t: 1* - (*p - t r )] t P -t r <t^t p 



Itr.sh. (t) 



(2) 



|tt [t - (t p - t d )] t p <t<t p +t d 



Then, the SF properties were studied with respect to the 
aforementioned time-scales and the main conclusion was 
that only the shortest time-scale between t r and td deter- 
mines the position of the break. Moreover, due to the fact 
that the ASCA data set of Mrk 501 consists mainly of sym- 
metric flares, they connect the shortest time-scale derived 
from the SF to the light crossing time of the emission re- 
gion. 

The apparent similarities between the observed SFs and 
those coming from the shot-noise simulations, presented by 
TAN01, give an erroneous impression that such correspon- 
dences hold. Apart from the fact that breaks in the SF can 
occur in data sets with no characteristic time-scales (Section 
[3j , we posit that the shot-noise model is not a physically re- 
alistic approach that can be used in order to associate the 
observed breaks of the SF to the smallest time-scales embed- 
ded in the data sets. In order for the shot-noise simulations 
to be realistic representations of the observed ASCA data 
set of Mrk 501, they must be able to reproduce the observed 
PSD as well as the SF. 

We repeat the simulations of TAN01 by creating shot- 
noise light curves 2000 t.u. long and we calculate the cor- 
responding PSD. Since the general shape of the latter can 
be understood qualitatively by the PSD form of the individ- 



ual triangular shots, used for the construction of each light 
curve, it is useful to look at the general analytical form of the 
PSD for a single triangular shot as derived directly from the 
Fourier transform of the continuous functions (Appendix [C)) . 
We should note here that since the simulated light curves 
are discretised, apart from their noisy form, due to the ran- 
domness of the occurrence times and of the intensities, we 
expect the PSD to flatten towards high frequenc ies due to 
aliasing effects (e.g. iPapadakis fc Lawrencelll993r ). 

Initially TAN01 consider that the light curve consists of 
identical symmetric triangular shots i.e., r — t r = td = W 
t.u. The PSD of such a process breaks at r _1 and becomes 
very steep at higher frequencies with power-law of index 
a w — - 4 superimposed on a wave with peaks separated by 
r _1 (Fig. [6] panel a). Qualitatively we can justify this be- 
haviour by considering the PSD of a single symmetric trian- 
gular shot which is proportional to /~ 4 sin 4 (7rr/) (equation 
(|C3p ). We expect that in the case of randomly occurring 
symmetric triangular shots, having different intensity but 
the same duration, the overall shape of their PSD should be 
similar to the sum of the PSDs of the individual symmetric 
triangular shots. 

The second case considered by TAN01 is that of non- 
identical symmetric shots (r = t T = td) whose time-scale r 
varies randomly between r m i n = 10 t.u. and r max = 100 t.u. 
In this case the PSD is flat until r~ 4 x an d then it breaks 
following a simple steep power-law of index a « —4 (Fig. 
[S] panel b). The oscillatory behaviour, caused by the sine 
term of equation (|C3|I . is smoothed out due to the random 
distribution of rs coming from every symmetric shot. 

The third and the last case considered by TAN01 is that 
of asymmetric identical shots with r r = 10 t.u. and Td = 100 
t.u. This is actually a generalization of the first case. The 
PSD breaks at t^ 1 , then it becomes steep with a power-law 
index of a w — 2 and after r r _1 beats separated by r r are 
superimposed on the underlying power-law (Fig. [5] panel c). 
Qualitatively, based on the PSD of a single asymmetric shot 
(Fig. IC1|I we can readily distinguish very similar spectral 
features such as the position of the first break, the slope 
between t^ 1 and r" 1 , and the separations in frequency of 
the beating frequencies. 

Each of the aforementioned simulated PSDs differs sig- 
nificantly from the PSD derived from the ASCA observa- 
tions i.e. a simple power-law with index a — —1.80 ± 0.09. 
Neither steep slopes nor sinusoidal features in oscillating or 
beating forms appear in the featureless PSD of Mrk 501. 
The fact that the PSD functions differ significantly between 
simulations and observations clarifies in an unambiguous 
way that the variability properties of Mrk 501 can not be 
described in a physically correct way from the shot-noise 
model. The latter can very well reproduce the various SF- 
features i.e. slopes and breaks observed in the actual SFs of 
the observed in blazars, but it can not reproduce the feature- 
less PSD. This means that the shot-noise simulations are not 
representative of the true underlying variability properties 
of Mrk 501 and hence should not be used in order to as- 
sociate the observed SF-breaks with characteristic physical 
time- scales. 
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Figure 6. Example PSDs for the shot-noise model consisting of 
triangular shots for three cases of shots. 

[Top panel (a)] Identical symmetric shots with rise and decay 
time, T r and Ty respectively, equal to 10 t.u. 

[Middle panel (b)] Symmetric shots with rise and decay time equal 
to t (tv = Td = t), and r distributed between r m i n = 10 t.u. and 

Tmax = 100 t.U. 

[Bottom panel (c)] Asymmetric identical shots with rise and decay 
time r r = 10 t.u. and Td = 100 t.u., respectively. 



5 FITTING MODELS TO THE STRUCTURE 
FUNCTION: PARAMETER AND ERROR 
ESTIMATION 

5.1 Lack of Gaussianity and statistical 
independence 

One of the major problems that affects the use of SF is 
that the various estimates, SF(ri), are not independent of 
each other. This problem severely affects the fitting routines 
e.g. least-squares, maximum likelihood, that are commonly 
used in the published blazar-SF-literature to derive the SF- 
breaks and the -slopes. These routines all require that the 
data points should be statistically independent so that the 
probability of obtaining the full ensemble, in our case SF(ri) 
for all Ti , is equal to the product of the probabilities of ob- 
tainin g the individual SF(ri) (e.g. iBevington fc Robinsonl 
Il992t ). 

That means that if we want t o fit a broken-power-law 
of the form (e.g. IZhang et aflbOO^ 



SF m (t; C, r max ,/3i,/3 2 ) = 




to a data set of the form 

{ (n, SF(n)) , (t 2 , SF(r 2 )) ,...,{r N , SF(r N )) } 

we have to maximize the probability of obtaining the en- 
semble of estimates {SF(ti), SFfa), . . . , SF(tn)}, known 
as the joint probability 

N 

P(C, W,/3i,/32) = nP t {SF M {n;C,T n ^,l3 1 ,l3 2 )) (4) 

i=l 

where Pi(j3Fy!_(Ti; C> Ttnax, P2)) is the probability of ob- 
taining the estimate SF(Ti). Equation Q is valid (i.e. the 
joint probability equals to the product of individual proba- 
bilities) iff the estimates {SF(ti), SFfa), ■ ■ ■ , SF(tn)} are 
independent of each other. 

Unfortunately the adjacent SF-estimates are far from 
independent. In Fig. [7]we show the autocorrelation function 
(ACF) of the first 25 logarithmic SF-estimates (out of 137) 
of Mrk 501 (Fig. [TJ after detrending them by subtracting 
a fifth degree polynomial. Note here that this detrending 
is crucial since we want to check the linear relations be- 
tween the SF-estimates (e.g. short-term dependence), and 
not long-term linear-dependence which is due to the general 
shape of the SF. In Fig. [7] we quantify the degree of linear 
dependence of the SF-estimates for the slope, the break and 
the plateau regions. It is impossible to have linearly cor- 
related values which are independent^. We can readily see 
that up to r = 0.88 days (the range which is used in order 
to derive the break time-scale in Section I3.1|l the degree of 
linear correlation is greater than 0.5. 

The assumption of Gaussianity is another issue with 

5 We should note that sometimes the absence of a linear correla- 
tion is confused with statistical independence but this is the case 
only when we are dealing with Gaussian distributions. Genuine 
statistical independence requires not only linearly uncorrclatcd 
values but also the absence of any functional relation between 
the variables. 
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Figure 7. The ACF of the first 25 detrended logarithmic SF- 
estimates of Mrk 501 (from Fig. [TJ showing the degree of linear 
dependence among them. The vertical dotted line shows the po- 
sition of the SF-break at 0.88 days, and the arrows indicate the 
SF-slopc and plateau regions. The dashed lines among the points 
of the ACF are linear interpolations intended to guide the eye. 



the SF fitting procedures. Only if the distribution of the 
SF-estimates within each ti is Gaussian (with a mean value 
SF(ji) and standard deviation a t ), is it valid to consider 
that 



1 r SF(T i )-SF M (T i ;C,T max ,ff 1 ,/3 2 )l 



27T 



(5) 



and estimate the product in the joint probability (equation 
Q) (which implies independence) as 



P(C,T max ,Pl,P2) = 
N 



n 



i 



1 y^ jv 



•SF(T < )-fF fcl (rj;C.Tn,ax.gl.g2)1 



(6) 



According to the maximum likelihood method, in order to 
find the most probable set of parameters (C, r max , /3i, P2) we 
must maximize P(C, r max , /3i, P2) (equation ([6])) or equiva- 
lently we must minimize the sum-argument x 2 in the expo- 
nential of equation ([6]) 



X 



SF(n)-SF M (n;C 1 'max 5 



(7) 



Thus the quantity \ automatically defines the goodness-of- 
fit which is actually a measure of the probability of obtain- 
ing the observed SF-estimates from a given set of parameters 
(C, Tmax, Pi, P2) but only when we are dealing with indepen- 
dent and Gaussian SF-estimates. 

However, in reality we do not deal with Gaussian SF- 
estimates. By producing 2000 random data sets 500 t.u. long 
from a PSD of a simple power-law form having index -2, we 
form the distribution of the SFs for the first, the second, and 
the three- hundredth time bins respectively (Fig. [8] top pan- 
els). The histograms clearly have a non-Gaussian form and 



remains non- Gaussian for all the time bins r. That means 
that any fitting relation of the form of equation @, deal- 
ing directly with the SF-estimates, does not provide reliable 
results since Gaussianity i.e. equation ([5]), is not valid. 

Interestingly, at least for those time bins which are be- 
low the fake SF-break, occurring at Tb r ~ 50 t.u. (Fig.0, the 
logarithms of the SF-estimates follow a distribution which 
is closer to Gaussian (Fig. [8] bottom panels a and b). For 
t > Tbr even the logarithmic-estimates of the SF deviate 
significantly from Gaussianity (Fig. [8] bottom panel c) and 
of course for all the time bins r the measurements continue 
to b e statistically depend ent on each other. 

IKataoka et al.1 (|200ll ) define the sum of square differ- 
ences X lm = J2k{^g 10 [{SF(r k ))] - log 10 [SF(r fe )]} 2 , where 
the angle brackets indicate the arithmetic mean, in order 
to derive a statistical significance of the goodness of their 
SF fit, accounting in this way for the non-Gaussianity. They 
argue that Xsim ls different from the traditional \ 2 (i- e - equa- 
tion (7J) but the statistical meaning is the same. Since the 
SF measurements are not independent, a \ 2 function of the 
form J2{S-E) 2 , where S comes from simulations (i.e. from 
a given assumed underlying model) and E are the actual 
estimates, makes sense only if for every SF bin the distri- 
bution of the entries is Gaussian (yielding the (S — E) from 
the exponent of the Gaussian, equation ((5}) and if all the SF 
bins are independent of each other (yielding the sum from 
the joint probability, equation (©). 

Conc erning the Gaussianit y it would be more appro- 
priate for IKataoka et al.1 l|200ll ) to use the mean logarithm 
of the SF, {\og 10 [SF(T k )]) rather than the logarithm of the 
mean SF, \og 10 [{SF(Tk))], in their pseudo chi-square esti- 
mate: x s 2 im , Gauss = J2 k {{log 10 [SF(r k )]) - log 10 [SF(r k )]} 2 
where both quantities (\og 10 [SF(T k )]} and \og 10 [SF(r k )] are 
distributed Gaussian within a tk- However, the main prob- 
lem of non-independence is not avoided and therefore the 
derived significances are not statistically meaningful. 



5.2 Error underestimation 

In most of the cases, when we deal with statistically inde- 
pendent variables under the assumption of Gaussianity, the 
most probable value (i.e. the one that minimizes the \ 2 ) an d 
its standard deviation are enough to give a full picture of the 
distribution of the fit parameters. 

There are several methods, as we are going to see, that 
are used in the literature claiming robustness and statisti- 
cally meaningful errors for the SF fit parameters i.e. slopes 
and breaks. In this section we show through simulations 
that the lack of statistical independence in the SF-estimates 
yields very small estimates of uncertainty in the fitting pa- 
rameters, which do not reflect the true underlying distribu- 
tions of the fitting parameters. The essence of any correct 
fitting procedure is to yield a statistically meaningful de- 
scription of these distributions. 



5.2.1 Fitting the SF-slopes and -breaks 

We use the previously generated simulations (2000 light 
curves, 500 t.u. long from a power-law PSD having an in- 
dex -2) , and for each light curve we calculate the logarithm 
of the SF-estimates. We attribute to every estimate a stan- 
dard error based on the error of the sample mean for each 
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Figure 8. The distributions of the ensemble of 2000 SF-estimates, coming from 2000 random data sets 500 t.u. long having a PSD of a 
power-law form of index -2, for the first, second and three-hundredth time bin respectively. 

[Top panels] The distributions of the normal SF-estimates (the histogram bins have a length of 1,1 and 100 respectively) is not Gaussian 
for all the r. 

[Bottom panels] The distribution of the logarithmic SF-estimates (the histogram bins have a length of 0.025, 0.025 and 0.05 respectively) 
approximates adequately the Gaussian distribution only for r < 50 t.u. 



time bin, which is one of the most usual methods (e.g. 
Collier fc Peterson! l200ll ; IZhang et al] 12001 ICzernv et all 
20051 ). We then fit, using a least-squares fitting procedure, 
a simple linear model of the form 



log 



; 10 [SF(r;C J /3)] = C + y 81og 10 (r) 



(8) 



where (3 represents the slope of the SF. Hence, for each light 
curve we estimate the value of the SF-slope (3 together with 
an error 8/3. The top panels of Fig. [9] show the distribution 
of f3 and 8/3. The errors derived from the fit, 8/3, are very 
small in comparison to the actual scatter of the measured (3, 
depicting clearly the effect of non-independency. From the 
simulations we expect 68 per cent of the derived slopes to 
be within the interval 1.02 ± 0.13 but erroneously the 8/3 
distribution has a mean of 0.006, dictating a range for the (3 
within 1.020±0.006. The reason for the very small uncertain- 
ties in the fitted slopes, 8(3, is that the errors on the actual 
SF-estimates are very small themselves due to their non- 
independent nature. Every SF-estimate follows smoothly the 
increasing trend, defined from the adjacent points, and thus 
the errors derived directly from the SF-estimates in this way 
are always too small. 

On the contrary the binned logarithmic periodogram- 
estimates together with their corresponding errors do pro- 
vide us with statistically correct information about the be- 
haviour of the slope. For each one of the previous simu- 
lated light curves we fit the binned logarithmic periodogram- 
estimates with a simple linear model and once again we de- 
rive the slope a and the error 8a. We can see from the bot- 
tom panels of Fig. [9l that 68 per cent of the a-estimates are 
expected to fall within —1.98 ±0.18, in accordance with the 
distribution of 8a (having a mean of 0.17) coming directly 
from the fits to the binned logarithmic periodograms of the 
simulated light curves. A meaningful error analysis originat- 



ing from a fitting procedure, is correct only when its predic- 
tive character (i.e. in our case for the PSD: 68 per cent of 
the measurements should be within the range —1.98 ±0.18) 
coincides with the actual fluctuational outcome of the pro- 
cess (i.e. —1.98 ± 0.17). These simulations show us clearly 
that this is not the case for the SF errors, which are much 
smaller than the true spread of the SF-estimates for the 
same variability process. 

Occasionally au thors attemp t to ta ke account some of 
the SF problems. In lAgudo et al.l |2006l ). for the case of the 
intra-day variable blazar S5 0716±71 observed at 86 GHz, 
the authors use the interpolated version of the SF which 
takes into account only the e rrors caused by the interpola- 
tion (|Quirrenbach et al] [2000). By repeating this procedure 
for our simulations we get a mean error in the slope of 0.011 
which is very small compared to the expected 0.13. Another 
weakness of this method is that several points in the SF have 
zero uncertainty (i.e. the intersection point o f the forward 
and backward interpolat ions, see figure 4 in Agudo et al.l 
2006 and figures 24-40 in lQuirrenbach et al.ll2000l ) therefore 
they do not really have any statistical meaning. 

In iFuhrmarm et all <|2008l ), again for the case of 
S5 0716+71 but observed at radio frequencies, the authors 
make use of the term "variability timescale" which they re- 
late directly to the SF-break. To estimate its uncertainty 
they use the saturation level po and the two fitted parame- 
ters of the SF function a, /3, SF(t ) = ar^ , (see equation s 
(B.2) and (B.3) in the appendix of iFuhrmann et all [20081 ). 
We use the simulations of Section 13.11 in order to see how 
well we can specify the position of these spurious SF-breaks 
by applying the same methodology, despite the fact that the 
simulated light curves do not contain any real characteris- 
tic time-scale. The distribution of the spurious break time- 
scales is shown in the left panel of Fig. [2] (grey area) , having 



Structure function: caveats and problems 11 



140 

« 120 

o 

| 100 

Z 80 
o 

b 60 
I 40 
Z 20 




150 



100 



50 



<y8>=1.02 
0-^=0.13 



0.6 


0.8 


1.0 


1.2 


1.4 













<«>= 

(Tir 



-1.98 
=0.18 



-2.3 -2.2 -2.1 -2.0 -1.9 -1.; 



-1.7 -1.6 



150 



100 



50 



OLi 
0.00 



(5/3) =0.006 



0.01 



0.02 



0.03 



0.04 



0.05 



60 



140 
120 

r. 
U 

a. 100 
i 80 



60 
40 
20 



(Sa)=0.ll 



0.10 



0.15 



0.20 

6a 



0.25 



Figure 9. The distributions of SF- and PSD-slopcs estimated from an ensemble of 2000 artificial light curves having an underlying PSD 
of a power-law form with slope -2. 

[Top panels] (Left) The distribution of the SF slopes has a mean of 1.02 and a standard deviation 0.13 (the histogram bins have a 
length of 0.025). (Right) The distribution of the errors coming from the fit 8/3, having a mean of 0.006 (the histogram bins have a length 
of 0.001). 

[Bottom panels] (Left) The distribution of the PSD slopes a has a mean of -1.98 and a standard deviation 0.18 (the histogram bins have 
a length of 0.02). (Right) The distribution of the errors coming from the fit, 8a, having a mean of 0.17 (the histogram bins have a length 
of 0.005). 



a mean value of 0.92 days and a standard deviation of .09 
days. By applying the method of iFuhrmann et all (|2008l ) in 
the simulations of Section 13.11 we find a mean error for the 
breaks of 0.006 days which differs significantly from the ex- 
pected 0.09 days. Obviously with such small (and incorrect) 
uncertainties every feature in the SF can be considered as a 
significant time signature of a speci al source property, as it 
seen in the right panel of figure 7 in lFuhrmann et all (|2008l ) 
where the authors find two "variability time-scales" embed- 
ded in their data set. 

Each of the above cases, show us that timescales de- 
rived from the commonly-used methods to interpret and fit 
the SFs in the blazar literature, should be treated with great 
caution. Due to the lack of statistical independence, the er- 
rors on the fitting parameters tend to be always very small 
underestimating their true statistical scatter. 

5.3 Resolving physically meaningful time-scales 
with the SF. 

True characteristic time-scales embedded in the blazar light 
curves, should be mapped in the SF as well as in the PSD. 
In the limit of a single PSD slope in the range -1 to -3, 
covering a large frequency range and for a long light curve, 
there is a direct theoretical relation between the PSD and 
SF slopes (see Appendix[B}. However for the more common 
broken power-law PSD and with typical astronomical ob- 
serving constraints, the relationships are less clear and thus, 
any attempt to relate SF-breaks or -slopes to PSD-breaks 



or -slopes by one-to-one relationships can lead to miscon- 
ceptions about the true underlying variability process. 

We produce 2000 artificial light curves 500 t.u. long, 
having a PSD of a broken-power-law shape with a break at 
/br = 0.1 (t.u.) -1 and with low and high frequency slopes 
of ct\ = — 1 and ah = —2.5 respectively. Then, for each light 
curve we estimate the SF and we produce an interpolated 
version of it in order to localize any feature around 10 t.u. in 
the form of either a local extremum or inflection point. We 
see a wide variety, and sometimes even an absence, of fea- 
tures around 10 t.u., something which places automatically 
the SF in the category of weak probing methods (Fig. Ill[) . 
In this particular case we ensure that any feature around 10 
t.u. is induced by the existence of the characteristic time- 
scale since we know a priori the expected position of it and 
hence no spurious break is expected at that time-scale. Of 
course this sort of interesting temporal signature has noth- 
ing to do with the previously mentioned (Section [3)l phys- 
ically uninteresting SF-break. As we can see from Fig. QT] 
these meaningless breaks, defined by the onset of a plateau 
on longer time-scales, are also present in our simulated light 
curves and their positions are simply deduced, as we showed 
in Section [3] from the length of the data set and the un- 
derlying PSD of the variability process. The distribution of 
the physically interesting SF-features, which are expected 
around 10 t.u., is shown in Fig. 1101 having a mean value of 
9.86 t.u. and a standard deviation of 1.63 t.u. coming from 
a total of 1789 SFs that exhibit features. 

In order to resolve these physically interesting breaks, 
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Figure 10. The distribution of SF-features (i.e. local extremum 
or inflection point) around 10 t.u., coming from 2000 artificial 
light curves 500 t.u. long, having a PSD of a broken-power-law 
shape with a break at /b r = 0.1 (t.u.)" 1 and with low and high 
frequency slopes of aj = —1 and a h = —2.5 respectively. The 
distribution has a mean value of 9.86 and a standard deviation of 
1.63. 



we employ the commonly used model of a broken-power-law, 
similar to the form of equation ([3]), but in its logarithmic 
version in order to ensure Gaussianity at least for the SF- 
estimates of the slope (Section 15. ip . 



loi 



ho[SF(r)] = 

C + Pi log 10 (r) T < T br 

C+(/3i- /3 2 ) log 10 (Tbr) + #a log 10 (r) t > r br 



(9) 



Fitting equation to the simulated data typically results 
in the minimisation routine finding the physically uninter- 
esting SF-breaks, tbr, occurring on long time-scales (~ 100 
t.u.) since this is the most prominent break-feature in the 
SF. Fortunately, and only for the case of the simulations, 
we know in advance the position of the expected SF feature 
we constrain our SF fits to the first 50 t.u., corresponding 
only to the increasing part of the SF, and we estimate the 
standard errors based on the error of the sample mean as 
described in Section 15.21 Moreover we constrain the range 
of rb r between 2-20 t.u. in order to force the fitted Tb r pa- 
rameter to correspond to the true PSD-break (Fig. 1111 solid 
line). 

The distribution of Tb r from the 2000 simulations, has a 
mean value of 5.31 t.u. and a standard deviation of 1.46 t.u. 
The shift of the most probable value from 9.86 t.u. to 5.31 
t.u. is clearly a systematic problem of that method which 
depends on the form of the fitted model (equation (0)). Nev- 
ertheless, the scatter of these breaks, 1.46 t.u., is very close 
to the expected one, 1.63 t.u. Once again, the errors on Tb r 
derived directly from the fitting procedure are extremely 
small, having a mean value of 0.40 t.u. The broken-power- 
law model is not very representative of the actual shape of 
the SF which is continuously curved. 



In order to model the physically interesting SF features 
we have to perform a much more detailed fitting even in 
this simple case tested here where we know in advance their 
approximate position. For this reason we em ploy a smoothly 
"dou ble-bending" power-law (section 4.2 in lMcHardv et all 
l2004h of the form 



SF(t) = Ct k 



1 + 



1 + 



T br- 



T 
Tbr, 



42 + ^1 



(10) 



with breaks at Tb ri and Tbr 2 an d power-law indices Ki for 
t < Tbrj, K2 for r bri < t < T br2 and K3 for r > T br2 . Again, 
we constrain the fitted break time-scales to be between 2-20 
t.u (Fig. [11] dashed line). 

For the 2000 simulated light curves the distribution of 
the fitted position of the first break, Tb ri , has a mean value 
of 6.13 t.u. and a standard deviation of 1.32 t.u., and for the 
second break 15.03 t.u. and 1.45 t.u. respectively. Neither of 
these breaks depicts the actual position of the characteristic 
time-scale (10 t.u.). Similarly to the simple broken-power- 
law fitting model, the errors derived directly from the fit are 
very small, having a mean value of 0.35 and 0.39 for the two 
breaks respectively. 

Despite the fact that breaks in the PSD, representing 
physically interesting variability time-scales, are mapped to 
the SF in the form of local maxima or inflection points, it 
is dangerous to make any statistical statement about them 
since 

(i) it is very difficult to depict the breaks robustly based 
on a fitting procedure due to the awkward shape of the SF; 
even if we know their position, based on the PSD, in which 
case we do not need to use the SF in the first place. 

(ii) the estimated errors in the fitted parameters (derived 
from a A\ 2 = 1, under the assumption of Gaussianity) 
are always very small with respect to the true scatter of 
the parameters due to the statistical dependence of the SF- 
estimates; even if we find a well behaved function able to 
depict the break based on single or multiple fitted parame- 
ters . 

In reality, we deal with astronomical data sets for which we 
want to reveal the possible existence of characteristic time- 
scales in them. In the case that these time signatures do 
exist, we do not know in advance their position and how they 
are mapped on the SF. Moreover, since the SF-estimates are 
affected by the statistical properties of the observed light 
curves (i.e. mean, variance and the length) and not from 
the actual statistical properties of the underlying variability 
process (as is the PSD), even in the absence of characteristic 
time-scales we expect to see SF-artefact-features similar to 
the ones induced by real characteristic time-scales. 



6 THE STRUCTURE FUNCTION AND 
DATA-GAPS 

The SF method is c onside re d by several rese a rchers 



(e.g. Takahash i et al. 120001; ICoflier fc Petersonl 1200 ll : 



Kataoka et al. 200l], |2002| ; IZhang et ail [2002) to be the 



ideal method of studying the time properties of gappy 
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Figure 11. The points represent the SF-estimates for a randomly 
selected light curve, drawn from the ensemble of 2000 data sets 
500 t.u. long, having a PSD of a broken-power-law shape with 
/br = 0.1 (t.u.)^ 1 , ai = —1 and = —2.5 respectively. The 
thick-dotted vertical line indicates the position of the expected 
characteristic time-scale of 10 t.u. The specific SF has an inflec- 
tion point at 10.32 t.u. A physically meaningless break appears 
at r ~ 100 t.u. 

The solid line represents the SF model- fit (equation J9}) em- 
ployed to depict Tb r which is constrained between 2-20 t.u., and 
yields Tt, r = 5.15 ± 0.43 t.u. The dashed line represents the SF 
model-fit (equation l(T0jl) employed to depict rt, r which is con- 
strained between 2-20 t.u., yielding T\, x . = 5.97 ± 0.38 t.u. and 



T"br 



14.92 ±0.37 t.u. 



The inlay shows the behaviour of the two fitted models around 
the expected characteristic time-scale of 10.32 t.u. 



data sets as it is believed to be less distorted by gaps than 
frequency-domain methods. Here, we use simulations to 
determine whether this confidence in the robustness to data 
gaps of the SF method is really justified. 

Initially we produce a single artificial light curve 2000 
t.u. long from a featureless PSD with a power-law shape with 
index -1.5 and we compute its SF (Fig. 1121 panels la and 
lb, respectively). We again note the spurious break at 315 
t.u. (as expected from Fig.JS]) which does not correspond to 
any real time property of the underlying variability process. 

We then add to the simulated light curve three different 
types of data-gaps and we estimate the SF (Fig. 1121 left pan- 
els). The first sampling scheme represents the case of almost 
periodic data-gaps where 57 per cent of the data are missing 
in almost periodic epochs, separated by 50-100 time units 
(Fig. 1121 panels lb). The second sampling scheme consists 
of dense and sparse sampling where data are equally spaced 
for small time periods, and more sparsely sampled for the 
remaining ones, yielding in total a data set having 83 per 
cent less data (Fig. 1121 panels lc). The last case that we 
consider is that of purely sparsely sampled data where 92 
per cent of the data are missing (Fig. 1121 panels Id). 

The effects of the abovementioned data-gaps on the SF 
are shown in the right panels of Fig. 1121 To take into account 



the uncertainties that are introduced by these gaps for the 
ensemble of our simu lated light curves, w e use the bootstrap 
method, discussed bv lCzernv et al.l (120031 ). From each gappy 
light curve we select randomly, allowing repetition, 1000 sub- 
samples of 2000 t.u. long, and we drop the multiple entries. 
Then we estimate the SF and from the ensemble of 1000 
SFs we calculate for every time bin r a standard deviation, 
including additionally the standard error in the sample mean 
coming from the gappy light curve. 

We should emphasise that here we study the differences 
between the SFs coming from the gappy data sets (empty 
diamonds, right panels of Fig. I12p and those coming from 
the uninterrupted data set (grey points connected with the 
grey line, right panels of Fig. [T2] which are also displayed in 
panel (lb) of Fig. I12|l , That does not mean that the unin- 
terrupted SF is necessarily the "true", "original", or "par- 
ent" SF of the underlying variability process. This is only 
one SF, chosen randomly from an ensemble of light curves 
having a featureless PSD of a power-law form with index 
-1.5. Our intention is to check whether the observed SF is 
affected by the d ata-gaps since this is the method that sev- 
eral authors (e.g . Takah ashi et alj 2000l: Collier fc Peterson! 



l200ll ; iKataoka et al.l l200ll . |2002l ; IZhang et all 12002ft use to 
derive physically interesting results. A direct visual compar- 
ison of the SFs (Fig. I12[) reveals that actually gaps do affect 
the shape of the SF since they introduce wiggles and bends. 
The latter can create the wrong impression that the under- 
lying variability process is described by an underlying PSD 
of a broken-power-law form (as in Section 

At this point it is interesting to check whether or not 
the SF-errors estimated from the bootstrap method depict 
correctly the deviations between the SFs of continuously 
sampled and gappy light curves, SFconti and 5Fgappy re- 
spectively. We produce 2000 artificial light curves 2000 t.u. 
long, having a power-law PSD of slope -1.5, and we esti- 
mate for each one of them the SF COTA i. Then, we apply to 
the light curves the dense and sparse sampling (example 
shown in Fig |12l panel 3a) and we estimate for each one 
of them the <S-F ga pp y as well as an error for each time bin, 
errboot,i{r) (i.e. the standard deviation of the gappy-SF of 
the i th simulated light curve at the time bin r, SF^ ppy (r)) 
as it is derived from the abovementioned bootstrap method. 
For each SF-bin we derive the following quantity 



|Alog 10 [SF(r) 
errboot,i(r) 



2000 1 , rorni 

1 y, [lQgiotggj, 



,(r)]-log 10 [S^ onti (T) 



2000 



err hoot ,i(T) 



(11) 



As we can see from Fig |13l the values of 
Alog 10 [S'F(r)]|err lm 1 ot A differ significantly from unity, 
having a mean value of 2.79 ± 0.91. That means that the 
bootstrap method does not yield statistically meaningful 
errors that reflect the true deviations between the con- 
tinuous and the gappy SFs. In other words, data gaps 
introduce systematic deviations in the structure function 
which depend on the light curve realization and cannot be 
accounted for without extensive simulations. 
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12. The effects of data-gaps in the SF of a simulated light curve (panel la) 2000 t.u. long, having a PSD of a power- 
5. The left panels show the gappy patterns. The right panels show the corresponding gappy SF-estimatcs (black points) 
continuous SF-cstimates (grey points connected with grey line), as they are shown in panel lb. 



law with 
together 



7 CONCLUSIONS 



We have performed an extensive series of simulations in 
order to test the properties of the classical "running vari- 
ance" SF-method. Under this fully controlled environment 
our studies have focused on the apparent SF-breaks, the 
correspondence of these breaks to the time-scales involved 
in the shot-noise model, the statistical behaviour of the SF 
with respect to the commonly used fitting procedures and 



its response to data-gaps. Our results can be summarized 
clearly as follows: 

• Strong SF-breaks frequently occur in data sets lacking 
any sort of characteristic time-scales. The position of these 
physically uninteresting breaks depends on the length of the 
observations and the shape of the underlying PSD. 

• SFs derived from blazar observations resemble those 
coming from the shot-noise model i.e. superposition of tri- 
angular shots. However in the frequency domain the same 
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Figure 13. The mean difference between gappy and continuous 
SFs coming from an ensemble of 2000 artificial light curves, hav- 
ing the same underlying PSD of a power-law form with index 
-1.5. The gappy sampling scheme, applied to the continuous light 
curves, is shown in panel (3a) of Fig. 1121 The errors err\, oot {r) 
are estimated for each light curve based on the bootstrap method 
(described in the text). The inset depicts the distribution of the 
points in the r = 1581 time-bin having a mean value of 4.16 and 
the standard deviation (i.e 34.1% of the entries) above and below 
the mean is shown in the grey region. 



model does not show the same behaviour as the observations 
and thus is not physically realistic. 

• Non-independence and non-Gaussianity mean that it is 
not possible to derive a meaningful goodness-of-fit from nor- 
mal fitting procedures. We see, for example, that the derived 
uncertainties on the SF quantities i.e. positions of breaks and 
slopes, are always much less than the actual scatter of these 
variables during multiple realizations of the same variability 
process. 

• Data-gaps affect severely the SF-estimates in an unpre- 
dictable way, introducing systematic deviations. The boot- 
strap method can not yield statistically meaningful errors 
depicting the deviations between the gappy and the contin- 
uous SFs. 

We finally comment that for blazar variability present im- 
plementations merely determine the shape of the SF for the 
one realisation under consideration. That one realisation is, 
of course, just one of many possible realisations and so does 
not define precisely the properties of the true underlying 
variability process, which is really what we want to know. 
Current PSD analysis methods (e.g. lUttlev et al1 l2002) gen- 
erally do, however, try to determine the shape of the un- 
derlying PSD by means of simulation-based modelling that 
take into account aliasing and other sampling effects. We 
encourage the blazar community to explore the use of these 
techniques or develop similar methodologies in the SF-field. 
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APPENDIX A: FUNCTIONAL RELATION 
BETWEEN THE STRUCTURE FUNCTION 
AND THE AUTOCORRELATION FUNCTION 

Consider TV observations of a stationary real-valued time 
process x(t): {x(ti), xfo), . . . , x(tjv)}, then there is a direct 
relation between the SF and the ACF of the process. The 
SF is given in general by 



1 



([x(t)-x(t + r)] 2 ) 



(Al) 



The autocovariance function of the time process x(i) is given 
by 



V x , x {t) = ( [a(t) - x(t)] [x(t + t) - £c(t)] ) 
The variance S of the process x(t) is equal to 

JV 



S 2 = 



i J2 [*(*o - *(*)] 2 = ( [*(*) - 



(A2) 



(A3) 



note that in the denominator we have N instead of TV — 1 
(i.e. the biased estimator) due to the fact that theoretically 
x(t) is estimated directly from the parent distribution and 
not from the data set itself. 
The ACF is given by 

Vx,x(t) 



ACF(t) = 



s 2 



(A4) 



By expanding the terms of equation ()A2|) and taking into ac- 
count the fact that we are dealing with a stationary process 
i.e. x(t) = x(t + t) = (x(t) \ = const, we get 



V x , x (t) = (x(t)x(t : + r)) - x(t) 
Similarly for the variance 



S 2 = x(t) 2 - x(t) = V X , X (0) 
From equation (|A1[) 



(A5) 
(A6) 
(A7) 



SF(t) = x(t) 2 - 2 (x(t)x(t + t)) + x{t + t) 2 
By means of equation (|A6[) 

SF(t) = 2 [s 2 + x(t) 2 - (x(t)x(t + r))] (A8) 

Finally from equation (|A4[l and equation (|A5|) the SF for a 
stationary process is given by 



SF(t) = 2 \S 2 



= 2S 2 



V x , x (t)] 
V x>x (t) 



s 2 



= 2S 2 [1 - ACF(t)] (A9) 
After equation (|A9() the normalized SF (NSF) is defined as: 

NSF(r) = SF{T) 



S 2 
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2 [l - ACF(t) 



(AlO) 



As t —> oo, ACF(t) — > (i.e. there is no linear correla- 
tion between the various measurements), from equation (|A9| 
SF(t) ->■ 2S 2 and from equation (|A"T0l) NSF(t) ->■ 2. 



APPENDIX B: FUNCTIONAL RELATION 
BETWEEN THE STRUCTURE FUNCTION 
AND THE POWER SPECTRAL DENSITY 

Consider a zero mean stationary real time process x(t) with 
autocovariance function (equation (|A5[n 



V x , x (t) = (x(t)x(t + t)) = {x(t)x*(t + r)) 



(Bl) 



where the asterisk denotes complex conjugation. By rewrit- 
ing the terms inside the mean as a function of their Fourier 
transform, Y(f), for — oo < / < oo, equation (|B1[) reads 

V.,*(t) = (/ Y(f)e- 2 ^ ft df Y'{f)e 2 ^ {t+T) df) 

\J — oo J — oo / 



Y(f)Y*(f)e MfT df 



\Y(f) e 



2irifT 



df 



(B2) 



and due to stationarity, none of the terms is varying as a 
function of t (that is, ind ependency of time translations e.g. 
iBendat fc Piersoll l| 19861 )) 



Vx,x(t) 



f 

J — ( 



\Y(f)\ e 



2irifT 



df 



(B3) 



According to the first equality of equation (|A9|) and equation 
161) for — oo < / < oo 



SF(r) 



2 [V x , x (0) - V x , x (t)] 



Y(f) df- 



Y(f) e 



df 



(1 



Y(f) df 



(B4) 



Based on the Euler's formula, e tK = cos(fi:) + isin(fi;), and 
ignoring the imaginary part dealing only with the phases 



SF(t) = 2 



[1 



s(27r/r)] \Y(f)\ df 



(B5) 



In order to estimate the variability power contained in the 
frequency interval between / and / + df, we can omit the 
distinction between positive and neg ative frequencies and 
regard / as varying between and oo (|Press et al.ll 19921 ) . By 
considering the sum of the modulus-squared of the sinusoidal 
amplitu des Y(f) and Y(— f), we estimate the one-sided PSD 
of x(t) (|Press et al.lll992i ) as 

-P(f) = \Y(f)\ 2 + \Y(-f)\ 2 for < / < oo (B6) 
Since x(t) G R, 

\ Y (f)\ 2 = \ Y {~f)\ 2 (symmetry about the y- axis) (B7) 
and thus, 

V{f) = 2\Y(f)\ 2 for < / < oo (B8) 



From the symmetric property (equation (|B7|) ) and due to 
the fact that cos(/) = cos(— /) for ^ / < oo, equation 
(|B5|) reads 

r°° i i 2 

SF(t) = 4 / [1 - cos(27r/r)] \Y{f) (B9) 
Jo 1 

yielding from equation (|B8|) 

r-OG 

SF(t) = 2 [1 - cos(27r/r)] T(f)df for ^ / < oo (BIO) 



For a PSD of a power-law form V(f) = k/~ a , with 1 < A < 
3 and k a positive constant the last integration (equation 
(|B10[0 has an analytical solution 



SF(t) = -2 A K7r A ~ 1 r(l - 



A) sin 



(Bll) 



where T(x) is the (complete) Gamma function. That means 
that only when the following requirements 

(i) stationarity (also called weakly stationarity i.e. mean 
value and autocova riance function indepe ndent of time 
translations, see e.g. IBendat fc Piersollll986h . 

(ii) zero mean data set. 

(iii) The frequency range / should vary from to oo. 

(iv) The PSD should be given from a power-law form with 
index 1 < A < 3. 

are fulfilled then there is a direct relation between SF and 
PSD connecting the slopes of the two quantities, /3 and A 
respectively, p = A — 1. 



APPENDIX C: THE GENERAL POWER 
SPECTRAL DENSITY OF A SINGLE 
TRIANGULAR SHOT 

The general PSD function of a triangular shot of the form 
of equation @ is given by 



-Px(f) = 



1-2 

to 



j , , j-4 {'■d + tdtr + t r 



87vH 2 t?f* 

(id + t T ) [tr cos(27rf d /) + id cos(27rt r /)] + 

t A t T cos(2n (t r + t d )f)} fori r / &t d + (CI) 

Figure [Cl] shows the form of this PSD function for the case 
of r r = 0.5 t.u., Td = 100 t.u. and Io = 0.04 flux units. 
It is characterized by two distinct breaks at f 1 — and 
fi — t t ~ and the slope changes from to -2 as we pass from 
/ < fi to f\ < / < /2 and becomes even steeper for f > fi 
with a slope of -4. For the latter region beat frequencies 
appear which are separated by A/ = r t . 

For instantaneous rise t T — the PSD is equal to 



-Px(f) 



;2 

J-0 



8 ^ 4f u 2/4 [i + ^ nf - coB(2irtd/) 

27rtd/sin(27rt d /)] 



(C2) 



for instantaneous decay the PSD is given again from the last 
relation (equation ()C2[) ') having the td replaced with i r . 

Finally, for the case of symmetric shots t T = t^ = r and 
equation (IC1|) reads 



Vx{f) 



[2 



7r 4 r 2 /' 



■ sin 4 (7rr/) 



(C3) 
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Figure CI. The general PSD function of a triangular-shaped 
shot (equation (in logarithmic scale). The shot parameters 
are r r = 0.5 t.u., = 100 t.u. and Iq = 0.04 flux units. The 
shape of the PSD is characterized by two breaks at /i = tJ -1 
and /2 = T r — 1 . For / < fx it is flat and for /i < / < /a it 
becomes steep with a slope of -2. For / > f% the PSD becomes 
even steeper with a slope of -4 and beat frequencies separated by 



A/i = 
trend. 



A/2 



A/n = T r start to appear in the general 



